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Abstract 


kh  e  convergence  rate  of  the  Schwarz  Alternating  Method  (fyM)  is  studied 
for  applications  involving  the  solution  of  elliptic  equations  on  composite 
grids.  Such  problems  arise  when  solvers  which  can  only  be  used  on  special 
domains,  such  as  rectangles,  are  used  for  more  general  region;  and  in  the 
disection  of  problems  for  parallel  processing.  It  is  shown  that  the  conver¬ 
gence  rate  is  a  function  of  the  overlap,  number  and  shape  of  the  subregions 
into  which  the  problem  domain  is  divided.  The  convergence  rates  for  fyM 
are  slow  and  an  accelerated  method  based  on  overrelaxation  techniques  is 
developed.  The  analysis  is  extended  to  predict  the  performance  of  the 
accelerated  method  and  optimal  relaxation  parameters.  Finally,  we  study 
the  effects  of  changing  the  iteration  order  for  the  fyM  and  accelerated  fyM 
methods 


<f,  -*v  *  ,  •  •  *  ,  •  «  •  •  •  «  '  O  \  *  */\*  »  *  -  '  •-* 


1  Introduction 


The  Schwarz  Alternating  Method  was  initially  introduced  to  obtain 

theoretical  results  for  elliptic  partial  differential  equations  [Sch69] .  This 
same  idea  has  subsequently  been  used  to  develop  algorithms  for  the  numer¬ 
ical  solution  of  the  same  class  of  problems  (see  [Tan86]  ).  This  work  has 
been  largely  limited  to  problems  in  two  dimensions  and  to  breaking  up  the 
problem  domain  into  limited  numbers  of  components.  The  mean  motiva¬ 
tion  for  using  the  method  has  been  to  deal  with  geometrical  difficulties — to 
reduce  the  problem  of  solving  on  a  complicated  domain  to  a  sequence  of 
problems  on  simple  domains.  This  has  often  been  motivated  by  a  desire  to 
use  special  solvers  that  can  only  be  used  on  regions  of  special  shape,  such  as 
rectangles.  These  results  confirm  the  applicability  of  the  method,  but  show 
that  the  original  iteration  is  limited  in  its  effectiveness,  its  convergence  rate 
is  slow. 

Our  interest  in  the  method  has  come  from  a  need  to  solve  elliptic  equa¬ 
tions  on  composite  grids  generated  to  both  resolve  complicated  geometries 
and  to  adaptively  resolve  the  solution  [CF085].  A  technique  is  needed  in 
this  situation  which  allows  one  to  use  a  standard  software  interface  for  each 
of  the  regions,  such  as  solving  on  a  sequence  of  rectangles,  and  which  min¬ 
imizes  the  communication  required  between  the  solution  processes  on  the 
subregions.  In  this  application  we  may  have  large  numbers  of  subgrids  in 
our  composite  grids.  We  are  also  interested  in  using  the  same  technique 
to  break  up  the  global  solution  process  into  loosely  coupled  subtasks  for 
parallel  processing.  This  can  also  lead  to  large  numbers  of  subgrids. 

We  have  begun  with  a  study  of  the  convergence  rate  of  the  original 
Schwarz  iteration  as  a  function  of  the  amount  of  overlap  of  the  subgrids. 
This  extends  the  earlier  work  of  Miller  [MiI65]  to  domains  in  Rp  and  to 
arbitrary  numbers  of  subgrids.  This  is  discussed  in  Section  2.  The  analysis 
is  carried  out  by  first  showing  that  Schwarz  iteration  can  be  formulated 
as  a  block  Gauss-Seidel  iteration  of  a  modified  system  of  equations.  The 
convergence  rate  is  found  to  be  exponential  in  both  the  amount  of  overlap 
and  in  the  number  of  regions. 

Since  the  convergence  rates  of  the  original  iteration  are  slow  and  lead  to 
inefficient  algorithms,  we  consider  a  modification  of  the  algorithms  where  we 
use  an  overrelaxation  technique  for  the  boundary  values.  This  iteration  is 
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described  in  section  3,  and  the  results  of  a  number  of  numerical  experiments 
are  presented.  It  is  found  that  the  strong  dependence  of  the  convergence 
rate  on  both  the  amount  of  overlap  and  on  number  of  subregions  can  be 
greatly  reduced  resulting  in  acceptable  performance.  It  is  demonstrated 
that  the  optimal  overrelaxation  parameter  cam  be  successfully  estimated 
using  a  formula  due  to  Young  [YouTlj.  The  sensitivity  of  the  optimal 
overrelaxation  parameter  as  a  function  of  overlap  and  the  number  of  grid 
parameters  is  also  studied.  Finally,  we  also  study  the  effect  of  changing 
the  order  in  which  the  iterations  are  carried  out  and  introduce  a  red-black 
iteration  which  can  be  used  for  parallel  processing  whereas  the  sequential  or 
Gauss-Seidel  iteration  cannot.  Fortunately,  it  is  found  that  the  red-black 
iteration  can  be  used  successfully.  Tables  of  data  which  illustrate  optimal 
choices  of  overlap  which  minimize  operations  are  also  included. 


2  Convergence  as  a  Function  of  Overlap 

The  numerical  Schwarz  algorithm  is  essentially  the  same  as  the  block  Gauss- 
Seidel  method  for  a  modified  matrix  equation  which  has  the  same  solution 
as  the  original  finite  element  or  finite  difference  equations  of  the  elliptic 
partial  differential  equation.  The  relationship  between  the  convergence  of 
Schwarz  Alternating  Method  (or  and  the  area  of  overlap  has  been 
observed  previously.  Miller  [Mil65]  proved  a  result  for  the  case  of  two  over¬ 
lapping  rectangles,  see  also  G.  Staxius  [Sta77]  and  E.  Volkov  [V0I68].  Miller 
was  mainly  interested  in  solving  elliptic  equations  in  irregular  regions,  an 
analysis  for  situations  which  arise  in  applications  motivated  by  parallel  pro¬ 
cessing  and  composite  grids  has  not  been  carried  out.  Here  we  will  show 
how  overlapping  affects  the  convergence  of  the  fyM  for  model  problems  in 
the  p-dimensional  case.  Before  we  consider  the  more  general  cases,  it  is 
easier  to  illustrate  the  analysis  for  a  one-dimensional  problem.  The  model 
problem  we  consider  in  one  dimension  is 

y"(x)  =  /(*),  *€(0,1) 

1/(0)  =  a;  y(l)  =  (3. 


Figure  1:  One  dimensional  overlapping  grid. 


After  discretization  using  a  centered  second  order  method  the  resulting 
linear  system  is 

Tnx  =  b,  (1) 

where 

Tn  =  Tridiagonal  { 1,  -2,  l}nxn. 

The  SfiM  for  solving  this  problem  divides  the  region  into  k  overlapping 
subregions  f L  i  —  1,  •  •  • ,  k  as  shown  in  Figure  1.  (To  simplify  the  analysis 
we  assume  the  overlapping  pattern  is  uniform,  similar  conclusions  can  be 
deduced  for  more  general  problems.) 

Let  h  be  the  grid  size,  t  the  length  of  the  overlap  and  r/  the  length  of  ev¬ 
ery  subregion.  Then  n  =  /  =  £  and  m  =  ^.  The  circular  points  in  Figure 

1  are  new  boundaries  for  the  subregions.  A  natural  way  to  implement  fyM 
is  to  first  guess  some  “reasonable”  initial  values  for  the  artificial  boundaries 
and  then  to  solve  these  subproblems  separately.  Next,  use  the  solutions  of 
these  subproblems  to  update  the  values  on  the  artificial  boundaries  and 
proceed  iteratively  until  the  solutions  on  the  overlapping  regions  converge. 
If  we  solve  on  these  subregions  in  a  natural  order,  each  succeeding  subre¬ 
gion  takes  its  new  boundary  values  from  the  new  solution  of  the  previous 
subregion.  This  procedure  is  equivalent  to  applying  the  block  Gauss-Seidel 
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method  to  the  following  modified  matrix  problem: 


=  (Tm  ®  Ik  +  Em  ®  Lk  +  Fm  ®  I7k)x  =  6. 


The  corresponding  block  Gauss-Seidel  iteration  for  this  equation  is  as  fol¬ 
lows: 

(Em  ®  I*  +  Tm  ®  J*)x(k+l)  =  — (Fm  ®  +  6. 

The  quantities  above  are  defined  as: 

•  Em:  an  m  x  m  matrix  with  zero  elements  everywhere  except  for  1  in 
position  (l,m  —  /  +  1). 

•  Fm:  an  m  x  m  matrix  with  zero  elements  everywhere  except  1  in 
position  (m,  /). 

•  /*:  a  k  x  k  identity  matrix. 

•  Lk'-  a  k  X  k  matrix  with  zero  elements  everywhere  except  for  l's  on 
the  subdiagonal. 

•  Uk'-  a  k  x  k  matrix  with  zero  elements  everywhere  except  for  l's  on 
the  superdi agonal. 

It  can  be  shown  that  (1)  and  (2)  are  equivalent  (see  [TanS6]).  Therefore, 
the  convergence  analysis  of  the  is  reduced  to  calculating  the  eigenvalues 
of  the  block  Jacobi  matrix  J  =  \I~l  N  where 

M  =Tm®  h , 

N  =  Em®  Lk  +  Fm®  Uk- 


If  we  multiply  out  M  lN  then 

J  =  (Tm  0  Ik)~l(Em  ®  Lk  +  Fm®Uk) 

=  (T"1  0>  Ik)(Em  0  Lk  +  Fm®  Uk) 

=  {T-lEm)®Lk  +  (T-lFm)®Uk 
=  Em  0>  Lk  +  Fm  0  Uk, 

where  Em  and  Fm  have  almost  all  zero  elements  except  columns  (m  —  /  +  1) 
and  l ,  respectively.  It  is  clear  that  the  rank  of  this  matrix  is  at  most  2k. 
After  some  row  and  column  exchanges  J  can  be  transformed  to  J,  which 
is  similar  to  J: 


J  —  U  JUT  —  ®(>»-2fc)x(n-2A:)  C(n- 

02fcx(n-2Jfe)  ^2, 


2k)x2k 


where 


G  =  D'  ®  Ik  +  E'  (&  Lk  +  F'  ®  Uk, 

Ol  n,_[0  6  1  w  _  [  0  0  ' 

L  0  0  ’  U  [6  0  J  ’  ^  ~  L  0  a  J  ’ 

_  /  _  m  +  1  —  / 

(m  +  1)’  (m  +  1) 

It  is  cleax  that  A  j  €  (0  U  Ac),  where  A  j  and  Ac  are  the  eigenvalues  of 
matrices  J  and  G,  respectively. 

Theorem  2.1  If  a  <  0.5  then  A c  satisfies  the  following  equations: 

\q  +  2  *  a*  cos  9  *  Ac  —  a2  +  b2  =  0, 

A c  *  sin(A;  *  9)  =  a  *  sin ((k  —  1)  *  9) 
where  9  is  a  parameter  [Tan86j. 


p  =  max {|  Ac  |}  =  max{|  \j  |}, 
it  is  easy  to  show:  if  k  =  2 

P  =  b, 

if  k  =  3 

p  =  y/b. 

Now  we  can  immediately  observe  some  important  facts  about  : 


OVERLAP  PUUIOTKR  (l/m)  ss  1  -  6 


Figure  2:  Theoretical  and  computational  values  of  the  squared  spectral 
radius  for  the  block  Jocobi  iteration  matrix  in  the  1-D  case.  Increasing 
k  =  l/m  corresponds  to  increasing  overlap. 

1.  First,  the  spectral  radius  of  J  only  depends  on  the  number  of  subre¬ 
gions  k  and  the  overlapping  area  a.  This  means  that  the  convergence 
of  is  independent  of  mesh  size  h. 

2.  For  the  case  of  k  =  2  and  3,  we  notice  that  when  the  overlapping 
area  increases,  p  decreases,  and  when  k  increases  p  increases.  These 
conclusions  also  are  valid  for  the  general  case  (k  >  3).  We  cannot  give 
a  closed  form  solution  for  k  greater  than  5  but  the  numerical  results 
indicate  that  similar  results  hold(see  Figure  2). 


upwards  to  improve  visibility  of  the  overlapping  pattern 

A  similar  analysis  can  be  applied  to  the  model  problem  in  2  dimensions. 
The  Poisson  equation  in  two-dimensions  is: 


d2U  d2U  ,,  , 

d^  +  ~a^  =  f{x-y)’ 


(x,y)  €  (0,1)  x  (0,1). 


Using  central  differences  we  obtain  a  discretization  of  this  equation: 


where 


Ax  =  b, 


A  =  Tn  ®  In  +  In  ®  Tn. 


This  is  of  the  same  form  as  we  obtained  in  the  one-dimensional  case,  h  is  the 
mesh  size  and  n  =  1.  If  we  cover  (0,1)  on  the  x  axis  with  k  subregions  as  in 
the  one-dimensional  case,  then  the  solution  area  is  covered  by  k  overlapping 
rectangles  as  shown  above. 


If  we  apply  to  these  overlapping  subregions,  then  it  is  equivalent  to 
applying  the  Gauss-Seidel  method  to  the  following  modified  equation: 

W  F'  1 

rm 

E'  IV  F' 

•C-m  m  rm 

Ax  -  X  (5) 


=  1  V-  •:  L 


E'  \V  F' 
E'  IV 


)  0  I*  +  (I,  3 Fm)  0  Uk}x  =  b, 


where 

=  T.,  *  /..  7.  E'm=EZ  Em,  fl  =  U  Fm. 

In  order  to  analy/e  *he  .  .r.v»Tg»*nre  of  w  ,  we  need  to  study  the  spectral 
radius  of  the  bio.  k  Ja.-obi  .terative  matrix  derived  from  the  matrix  A: 


j  = 


where 


•V  -  11m  3  A  =  {In  3  Em  )  3  i)c  +  (In  0  Fm)  0  Uk- 

Then  we  have  the  following  result: 

Theorem  2.2  The  matrix  J  is  similar  to  the  matrix 


fin3 — 2nk X n2 — 2nk  G  nl ~2nkx2nk 


* 2nk  x  n*  —  2nk 


r2nk  X  2 nk 


G  =  Block  —  diagonal  {G,} ,  i  =  1,  ••  • ,  n  , 

G,  =  D[  0  Ik  +  E[  2>Lk  +  F:  0  Uk, 


a,  = 


cr  -  cr 


3,  = 


>  ( 1  — ^)m  —  (1  —  *)m 

St  Si 


cr  -  cr m  cr  -  crm 

C.  =  V,  +  \jv}  -  1,  Vi  =  2  -  cos(  ;7r  ), 
v  n  +  1 

and  k  —  l / m  is  the  ratio  of  the  overlap. 


i  =  1, • ■  •  ,  n 
i  =  1, • • •  ,  n 


Proof.  Let 


U  =  (An  S  Im)  ®  J* 


where  An  is  an  orthogonal  matrix  for  which  the  columns  are  the  eigenvec¬ 
tors  of  the  matrix  Tn,  and  U  is  an  orthogonal  matrix.  Note  UNUT  —  N . 


Then 


r  =ujut  =  um~1  nut 

=  (UMUt)-1N 

=  {((AnTnAj  ®  Im)  +  In®  Tm)-1  ®  Ik}N 
=  {{Dn  <g>  Im  +  In  ®  Tm)_1  ®  /fc}lV  , 

where  Dn  is  a  diagonal  matrix  whose  diagonal  elements  axe  the  eigenvalues 
of  T„.  We  know  that  there  is  an  ran  x  mn  permutation  matrix  P  such  that 

P{A®B)Pt  =  B®  A 

where  A  and  B  are  any  n  x  n  and  mxm  matrices,  respectively.  So  we  have 
P(/n  ®  Em)PT  =  Em.®  P(In  ®  Fm)PT  =  Fm  ®  In, 

P(T»n  ®  /m  +  I»  ®  Tm)PT  =  Im  ®  £>„  +  Tm  ®  Jn. 

Notice  that: 


Q  =  Im®Dn  +  Tm®  In 

=  Block-diagonal{T,  }nx„, 

T,  =  Trzdmgona/{  1,7,,  l}m*m, 

Z  7T 

* ,  =4  +  2  cos( - ),  z  --  1,  ■  •  ,  n. 

n  4-  1 
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Let  P  =  P  <S>  Ik,  Then  we  have 

J"  =PJ'PT  =  P{(Drl®Im+In®Tm)-l®Ik}PTPVpT 

=  {(P(Dn  ®  Im  +  In®  Tm)PT)~l  ®  Ik}{{Em  ®  In)  S  Lk  +  (  /v,  0  In  )  ■:  l\ 
=  ( Q~x  ®  Ik){{Em  ®In)®Lk  +  (Fm®  In)  ®  Uk) 

=  ( Em  ®  In)  ®  Lk  +  (Fm  ®>  In)  ®  Uk 

As  in  the  one-dimensional  analysis,  we  can  move  all  of  the  non-zero 
columns  to  the  last  columns  and  the  theorem  follows. 

Since  the  structures  of  these  diagonal  blocks  are  the  same  as  those  ana¬ 
lyzed  in  the  one-dimensional  case,  we  can  find  a  similarly  tight  estimate  of 
the  spectral  radius  of  J,  pj,  by  using  Theorem  1.  But  here  it  is  clear  that 

a,  4-  &  <  1,  a,  >  0,  >  0  i  =  1,  •  •  • ,  n 

and  thus  we  cannot  derive  a  closed  form  for  pj,  but  we  may  use  the 
Gershgorin  theorem  to  get  a  bound  for  pj. 


Corollary  1 


Pj  <  cti  +  0\- 


If  we  denote  p  =  it  is  easy  to  estimate  the  asymptotic  bound  for  pj{ as 
h~>  0): 


Corollary  2  If  k  =  2, 


If  k  >  2 


sinh((l  —  K)pn) 
sinh(/i7r) 


sinh(«^7r)  +  sinh((l  —  K)pn) 
sinh(^i7r) 


The  relationship  between  pj  and  the  two  parameters  /  and  p  is  shown  in 
Figure  4. 

Similar  results  can  be  derived  for  model  problems  in  more  dimensions. 


>  2  SUBREGIONS 

=  1.0 


/I  =1.0 

3  SUBREGIONS 


/i  =  2.0 


3  SUBREGIONS 


OVXRULP  PARAMETER  ( 


Figure  4:  Theoretical  bounds  and  computational  values  of  the  squared  spec¬ 
tral  radius  for  the  block  Jocobi  iteration  matrix  in  the  2-D  case.  Increasing 
i  corresponds  to  increasing  overlap.  Note  that  the  domain  size  increases 
with  increasing  overlap  when  the  subregion  size  is  held  fixed. 


Wrl  ‘1 
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Theorem  2.3  For  the  p-dimensional  model  problem  the  asymptotic  bound 
for  the  spectral  radius  of  the  block  Jacobi  iterative  matrix  is: 

—  l/c/iT r)  +  sinh(\/p  —  1(1  —  K)pr) 
sin  h(v//>  —  Ia471") 

Here  we  use  the  same  notation  as  used  for  the  one-dimensional  case:  l  =  ~ , 
m  =  ?,  k  =  —  and  u  =  —. 

The  same  analysis  can  be  applied  to  the  nine-point  stencil  or  tensor 
product  finite-element  schemes  which  have  higher  accuracy.  The  latter 
need  some  more  complicated  eigenvalue  analysis.  We  will  not  list  all  of 
these  results  here,  but  our  results  indicate  that  the  relationship  between  the 
convergence  of  and  the  overlapped  area,  and  the  relationship  between 
convergence  and  the  number  of  the  subregions  is  similar.  All  these  results 
guide  us  to  more  efficient  implementations  of  the 


^  sinh(  s/p 
PJ  <  - 


m 


3  Numerical  Results  and  Over  Relaxation 


Numerical  experiments  were  performed  to  verify  the  theoretical  results  of 
the  previous  section  and  also  to  search  for  methods  which  more  efficiently 
solve  the  overlapping  grid  problem.  The  1-D  model  problem  (1)  and  the  2- 
D  Poisson  equation  (4)  are  solved  numerically  using  the  SfiM  ■  Convergence 
of  the  numerical  solution  to  the  exact  solution  is  measured  by  calculating 
the  l\,  I?  and  /«,  norms  of  the  solution  error.  These  norms  are  defined  by 


l|tf  lit  =  II  1  “  u'-j  I  AxAy 


ll^lh  = 


J2  I  -  u'o  I2  AxAy 

ll^llco  =  max  I  U,J  -  u,tJ  I 


1/2 


where  is  the  exact  solution. 

We  solved  the  one  and  two  dimensional  problems  calculating  the  error 
norms  at  each  iteration  of  the  .  The  solution  error  norms  all  decay 
exponentially  at  the  same  rate,  this  decay  being  of  the  form 


\ew\\i  =  lk(O)ll.0fc 
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where  k  is  the  iteration  number.  The  constant  o  is  the  convergence  factor 
and  estimates  of  it  should  compare  with  the  theoretical  results  for  the 
spectral  radius  of  the  G-S  iterative  matrix  (2).  The  theoretical  estimates 
of  the  previous  section  were  for  the  iterative  Jacobi  matrix  but  the  G-S 
spectral  radius  is  the  square  of  the  Jacobi  spectral  radius  (Varga,  1962). 
Hence,  we  can  directly  compare  the  theoretical  and  experimental  results. 

Figure  3  is  a  plot  of  the  theoretical  estimates  of  the  spectral  radii  of  the 
G-S  iterative  matrix  and  the  experimental  convergence  factor  versus  the 
nondimensional  overlap  l/m  for  the  1-D  problem  with  2,  3,  4  and  6  sub- 
regions.  Theory  and  computation  compare  well,  especially  at  the  smaller 
overlaps  (larger  values  of  b).  As  the  theory  indicates,  the  Jacobi  spectral 
radii  provide  an  upper  bound  on  the  convergence  factors  and  the  conver¬ 
gence  factors  increase  with  smaller  overlaps. 

Two  dimensional  problems  are  characterized  by  two  parameters, 
H  and  l.  Corollary  2  gives  expressions  for  the  spectral  radii  of  the  two 
subregion  iterative  matrix,  and  also  an  upper  bound  on  the  spectral  radii 
for  more  than  two  subregions.  The  theoretical  prediction  of  the  spectral 
radii  for  the  G-S  iterative  matrix  along  with  numerically  determined  decay 
rates  for  the  2-D  problem  are  shown  in  Figure  4.  The  relative  grid 
ratio  n  is  held  constant  and  the  relative  overlap  /  is  varied  for  the  2,  3, 
4  and  6  subregion  SaM  problem.  Once  again  computations  and  theory 
agree.  The  behavior  in  the  2-D  case  is  similar  to  that  of  the  1-D  case  with 
the  convergence  factors  increasing  for  smaller  overlaps.  The  theory  also 
correctly  predicts  the  behavior  of  the  convergence  factor  for  different  grid 
shapes  as  can  be  seen  in  the  agreement  between  theory  and  observations 
for  the  n  equal  to  two  case. 

The  monotonic,  exponential  decay  of  the  error  norms  indicates  that 
some  form  of  overrelaxation  may  be  suitable  for  increasing  the  convergence 
rate  of  the  numerical  solution  to  the  true  solution.  One  possible  scheme  is 
to  overrelax  at  the  boundaries  of  the  subregions.  If  u  is  a  boundary  value 
of  region  two  and  lies  in  region  one  then  we  could  set  this  boundary  value 
using 

u(‘)  =  u(‘-1)+w(u(t)-u(*-I))  (6) 

where  k  is  the  SfiM  iteration  step  and  u  is  the  value  obtained  from  the  latest 
solution  on  region  one.  The  overrelaxation  parameter  is  u?  .  u  =  1  produces 


no  relaxation  and  typically  1  <  w  <  2  This  can  be  thought  of  as  using 
with  SOR  (Successive  Over  Relaxation). 

The  results  obtained  using  the  overrelaxed  are  impressive.  Figures 
5  and  6  show  the  relative  work  used  to  obtain  an  error  of  less  than  10-5  (as 
measured  by  the  convergence  norms)  using  no  relaxation  factor  and  also 
using  a  optimal  relaxation  factor  for  the  1-D  -  two,  three,  four  and  six 
subregion  cases.  As  the  overlap  becomes  smaller  the  gain  using  the  over¬ 
relaxed  becomes  greater.  However,  the  total  amount  of  work  necessary 
to  obtain  a  converged  solution  is  now  only  weakly  dependent  on  overlap. 
Table  1  lists  the  overrelaxation  parameters  used  in  these  cases.  Figures  7 
and  8  along  with  Table  2  contain  results  for  the  2-D  model  problem.  For 
the  2-D  problem,  all  the  points  along  an  inner  boundary  are  overrelaxed  us¬ 
ing  (6).  The  2-D  results  do  not  show  as  dramatic  savings  in  computational 
effort  as  are  found  in  the  1-D  case;  but  the  savings  are  still  very  significant, 
usually  cutting  the  amount  of  work  in  half. 

As  in  other  applications  of  relaxation  techniques,  the  critical  task  is 
choosing  the  relaxation  parameter  omega.  Figure  9  is  a  plot  of  the  relative 
work  versus  the  value  of  omega  for  various  1  and  2-D,  n  subregion 
problems.  As  the  amount  of  overlap  decreases  the  computational 
savings,  optimal  omega  and  the  sensitivity  of  the  savings  to  the 
choice  of  omega  all  increase.  The  sensitivity  of  the  computational 
savings  to  the  choice  of  omega  is  much  smaller  in  the  2-D  problem  compared 
to  the  1-D  problem.  As  we  would  expect,  the  choice  of  omega  is  most  critical 
when  the  possible  savings  are  greatest. 

Since  matrix  A  in  equation  (5)  has  the  property  A^[You71],  we  can 
use  the  following  formula  for  the  calculation  of  an  optimal  overrelaxation 
factor  given  the  spectral  radius  of  the  Jacobi  iterative  matrix. 

2 

w opt  = 

Also  in  tables  1  and  2  are  the  optimal  relaxation  factors  calculated  from  this 
formula  along  with  those  found  experimentally.  For  the  1-D  problem  the 
prediction  of  omega  is  very  good  and  fortunately  is  best  when  the  overlap 
is  smallest,  i.e.  when  there  is  the  most  to  gain  and  the  choice  of  relaxation 
factors  is  most  critical.  In  the  2-D  problem  the  prediction  of  omega  is  not  as 
good  as  in  the  1-D  problem,  but  the  sensitivity  of  the  savings  to  the  choice 
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Figure  5:  Theoretical  and  computational  estimates  of  the  relative  work 
for  the  1-D  case  using  no  relaxation.  The  work  is  relative  to  that  needed 
to  solve  once  on  a  single  grid  over  the  entire  domain.  Solution  work  is 
proportional  to  the  number  of  points  on  a  grid. 


of  omega  is  much  less  in  the  2-D  problem,  and  again,  as  in  the  1-D  case, 
the  prediction  of  the  optimal  relaxation  factor  is  better  as  the  sensitivity 
increases.  Almost  all  of  the  computational  savings  can  be  realized  using 
the  theoretically  predicted  optimal  relaxation  factors  in  both  the  1  and  2-D 
problems! 

Other  methods  of  increasing  the  convergence  rate  of  the  problem 
can  involve  changing  the  order  of  the  solution  of  the  subregions  for  cases 
where  there  are  more  than  two  subregions.  Two  possibilities  are  investi¬ 
gated.  The  first  involves  sweeping  forward  (as  in  SfiM  )  and  then  sweeping 
backwards  along  the  row  of  overlapping  grids,  and  hence  can  be  thought 
of  as  an  SSOR  (symmetric  SOR)  technique.  The  second  procedure  uses 
a  Black-Red'  solution  pattern,  solving  on  every  other  grid  during  the  first 
step  and  then  solving  on  the  remaining  grids  the  second  step. 

Sweeping  in  both  directions  (SSOR)  produced  little  to  no  gain  in  effi¬ 
ciency  compared  with  the  standard  ,  and  in  some  cases  requires  more 
computational  effort  to  produce  a  converged  solution.  When  combined  with 
overrelaxation  it  does  not  produce  as  dramatic  savings  as  does  the  with 
overrelaxation,  but  the  choice  of  relaxation  factor  is  less  critical.  There  is 
no  great  advantage  which  leads  one  to  choose  SSOR  over  SOR. 

The  Black-Red  solution  technique  produced  results  almost  identical  to 
the  results  for  both  the  unrelaxed  and  overrelaxed  methods.  The  com¬ 
putational  savings  and  optimal  relaxation  factor  are  the  same  for  a  given 
problem  and  hence  the  theoretical  determination  of  the  overrelaxation  fac¬ 
tor  can  be  used  with  the  B-R  scheme.  The  major  difference  between  the 
and  the  B-R  method  is  that  in  the  S$1  one  can  solve  on  only  one  subre¬ 
gion  at  any  time,  for  its  solution  is  necessary  to  set  boundary  conditions  for 
the  next  subregion.  Using  the  B-R  scheme  one  can  solve  on  up  to  half  the 
subregions  at  any  one  time,  each  one  being  completely  independent  of  the 
others  for  it  is  only  the  next  set  of  grids  that  needs  these  solutions  for  set¬ 
ting  boundary  conditions.  Thus  the  B-R  scheme  is  ideal  for  multi-processor 
machines. 
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discussions. 
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Figure  7:  Theoretical  estimates  of  the  relative  work  for  2-D  case  using  no 
relaxation  and  holding  the  domain  size  constant. 
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Figure  8:  Theoretical  estimates  of  the  relative  work  for  the  2-D  case  using 
optimal  overrelaxation  and  holding  the  domain  size  constant. 
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